#ifndef __NAVIERSTOKES_H
#define __NAVIERSTOKES_H 1

#include <life/options.hpp>
#include <life/lifecore/application.hpp>

#include <life/lifealg/backend.hpp>
#include <life/lifediscr/functionspace.hpp>
#include <life/lifepoly/im.hpp>
#include <life/lifefilters/gmsh.hpp>
#include <life/lifefilters/exporter.hpp>
#include <life/lifefilters/gmshtensorizeddomain.hpp>

#include <life/lifemesh/elements.hpp>
#include <life/lifevf/vf.hpp>
#include <life/lifediscr/bdf2.hpp>
#include <life/lifediscr/ale.hpp>



namespace Life
{

    template<int Dim,
             int Order,
             template<uint16_type,uint16_type,uint16_type> class Entity>
    class NavierStokes
    {
    public:

        typedef double value_type;

        typedef Application application_type;
        typedef boost::shared_ptr<application_type> application_ptrtype;

        typedef Backend<value_type> backend_type;
        typedef boost::shared_ptr<backend_type> backend_ptrtype;

        /*matrix*/
        typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
        typedef typename backend_type::vector_ptrtype vector_ptrtype;

        //___________________________________________________________________________________//

        /*mesh*/
        typedef Entity<Dim,1,Dim> convex_type;
        typedef Mesh<convex_type> mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

        //___________________________________________________________________________________//

        /*basis fluid*/
        typedef Lagrange<Order, Vectorial> basis_fluid_u_type;
        typedef Lagrange<Order-1, Scalar> basis_fluid_p_type;
        typedef bases<basis_fluid_u_type,basis_fluid_p_type> basis_fluid_type;

        /*space fluid*/
        typedef FunctionSpace<mesh_type, basis_fluid_type> space_fluid_type;
        typedef boost::shared_ptr<space_fluid_type> space_fluid_ptrtype;

        BOOST_MPL_ASSERT( ( boost::is_same<typename space_fluid_type::bases_list, basis_fluid_type> ) );
        BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_fluid_type::bases_list,mpl::int_<0> >::type, basis_fluid_u_type> ) );
        BOOST_MPL_ASSERT( ( boost::is_same<typename mpl::at<typename space_fluid_type::bases_list,mpl::int_<1> >::type, basis_fluid_p_type> ) );

        /* functions */
        typedef typename space_fluid_type::element_type element_fluid_type;
        typedef boost::shared_ptr<element_fluid_type> element_fluid_ptrtype;

        typedef typename element_fluid_type::template sub_element<0>::type element_fluid_0_type;
        typedef typename element_fluid_type::template sub_element<1>::type element_fluid_1_type;

        //___________________________________________________________________________________//

        /* time */
        typedef Bdf<space_fluid_type>  bdf_type;
        typedef boost::shared_ptr<bdf_type> bdf_ptrtype;

        //___________________________________________________________________________________//

        /* export */
        typedef Exporter<mesh_type> export_type;

        //___________________________________________________________________________________//

        /* ALE */
        typedef typename space_fluid_type::template sub_functionspace<0>::type::element_type space_fluid_velocity_type;
        typedef typename space_fluid_type::template sub_functionspace<0>::type space_fluid_velocity_ptrtype;
        typedef typename space_fluid_velocity_type::element_type element_fluid_velocity_type;
        typedef boost::shared_ptr<element_fluid_velocity_type> element_fluid_velocity_ptrtype;

        typedef Bdf< space_fluid_velocity_type > bdf_ale_velocity_type;
        typedef boost::shared_ptr<bdf_ale_velocity_type> bdf_ale_velocity_ptrtype;

        typedef ALE< convex_type > ale_map_type;
        typedef boost::shared_ptr< ale_map_type > ale_map_ptrtype;

        typedef typename ale_map_type::pN_functionspace_type ale_map_functionspace_type;
        typedef typename ale_map_type::pN_element_type ale_map_element_type;
        typedef boost::shared_ptr<ale_map_element_type> ale_map_element_ptrtype;

        typedef Bdf< ale_map_functionspace_type > bdf_ale_map_type;
        typedef boost::shared_ptr<bdf_ale_map_type> bdf_ale_map_ptrtype;

        // Interface
        typedef Mesh<Entity<1, 1,1> > mesh_structure_type;
        typedef boost::shared_ptr<mesh_structure_type> mesh_structure_ptrtype;

        typedef Lagrange< Order, Scalar/*, Continuous, double, PointSetWarpBlend*/ > basis_structure_type;
        typedef FunctionSpace<mesh_structure_type, basis_structure_type > structure_space_type;
        typedef boost::shared_ptr<structure_space_type> structure_space_ptrtype;
        typedef typename structure_space_type::element_type element_structure_type;
        typedef boost::shared_ptr<element_structure_type> element_structure_ptrtype;

        mesh_ptrtype createMesh( double meshSize );
        mesh_structure_ptrtype createStructMesh( double meshSize );

        NavierStokes( application_ptrtype __appli )
            :
            M_lifeApplication(__appli),
            M_backend( backend_type::build( this->application()->vm() ) ),
            meshSize( this->application()->vm()["hsize"].template as<double>() ),
            mesh_fluid(this->createMesh( meshSize ) ),
            mesh_reference_fluid(this->createMesh( meshSize ) ),
            mesh_struct(this->createStructMesh( meshSize ) ),
            Sh(new structure_space_type( mesh_struct )),
            // reference boundary description
            M_bc_reference_top(new element_structure_type( Sh, "bc_reference_top" ) ),
            M_bc_reference_bottom(new element_structure_type( Sh, "bc_reference_bottom" ) ),

            Xh(new space_fluid_type( mesh_fluid ) ),
            U_fluid( new element_fluid_type(Xh,"ub")),
            M_identity_ale(new element_fluid_velocity_type( Xh->template functionSpace<0>(), "identity_ale" ) ),
            M_mesh_velocity(new element_fluid_velocity_type( Xh->template functionSpace<0>(), "mesh_velocity" ) ),
            M_bdf_fluid ( new bdf_type( this->application()->vm(),Xh,"fluid" )),
            M_depl_top(new element_structure_type( Sh, "bc_top" )),
            M_depl_bottom(new element_structure_type( Sh, "bc_bottom" )),
            M_aleFactory( new ALE< convex_type >( std::make_pair(0,4), mesh_reference_fluid, this->application()->vm() ) ),
            M_aux(new ale_map_element_type( M_aleFactory->getDisplacement().functionSpace(), "new ale map" ) ),
            M_bdf_ale_map ( new bdf_ale_map_type(this->application()->vm(),
                                                 M_aleFactory->getDisplacement().functionSpace(),
                                                 "ale_map") ),
            M_bdf_ale_identity ( new bdf_ale_velocity_type( this->application()->vm(),
                                                            Xh->template functionSpace<0>(),
                                                            "velocity fluid" ) ),
            exporter( Exporter<mesh_type>::New( this->application()->vm(), this->application()->about().appName() ) )
         {
             M_cpt_export=0;
            //Masse volumique
             rho = this->application()->vm()["rho"].template as<double>();
             mu = this->application()->vm()["mu"].template as<value_type>();
             penalbc = this->application()->vm()["bccoeff"].template as<value_type>();
             Vitesse_Entree = this->application()->vm()["entree-v"].template as<value_type>();
             Pression_Sortie = this->application()->vm()["sortie-p"].template as<value_type>();

            *M_bc_reference_top = vf::project( Sh, elements(Sh->mesh()), vf::constant(1.0) );
            *M_bc_reference_bottom = vf::project( Sh, elements(Sh->mesh()), vf::constant(0.0) );
            M_referencePolyBoundarySet.push_back(*M_bc_reference_bottom);
            M_referencePolyBoundarySet.push_back(*M_bc_reference_top);

            // set up the non linear solver
            M_backend->nlSolver()->residual = boost::bind( &NavierStokes<Dim, Order, Entity>::updateResidual,
                                                           boost::ref( *this ), _1, _2 );
            M_backend->nlSolver()->jacobian = boost::bind( &NavierStokes<Dim, Order, Entity>::updateJacobian,
                                                           boost::ref( *this ), _1, _2 );
         }

        void run();
        void runPtFixe();

        void updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J /* , vector_ptrtype& R*/);
        void updateResidual( const vector_ptrtype& X, vector_ptrtype& R );
        void updateResidual2( const vector_ptrtype& X, vector_ptrtype& R );

        void updatePtFixe( element_fluid_type & U, element_fluid_type & Uold,element_fluid_velocity_type & mesh_velocity,
                           ale_map_element_type & depl );
        void AlgoNewton(vector_ptrtype & U);

        template<typename ExprType>
        void modifVec(element_fluid_velocity_type & u,ExprType expr);

        application_ptrtype application() { return M_lifeApplication; }

        void init();

        void updateDeplacement(element_structure_ptrtype __depl_top,
                               element_structure_ptrtype __depl_bottom );

        element_structure_ptrtype getDeplTop() { return M_depl_top; }
        element_structure_ptrtype getDeplBottom() { return M_depl_bottom; }
        void updateALEmap();

        void updateBdf();

        void exportResults() { exportResultsInstat(mesh_fluid,*U_fluid,*M_mesh_velocity,*M_aux, M_cpt_export); }

        void execNewton();
private:

        void solve( sparse_matrix_ptrtype const& D, element_fluid_type& u, vector_ptrtype const& F, bool is_sym );

        double diffL2( element_fluid_type & u1, element_fluid_type & u2 );

        void exportResults( element_fluid_type& u );

        void exportResultsInstat( element_fluid_type& u,
                                  element_fluid_velocity_type & mesh_velocity,
                                  ale_map_element_type & aux,
                                  double time );

        void exportResultsInstat( mesh_ptrtype __mesh,
                                  element_fluid_type& u,
                                  element_fluid_velocity_type & mesh_velocity,
                                  ale_map_element_type & aux,
                                  double time );


    private:

        application_ptrtype M_lifeApplication;

        backend_ptrtype M_backend;

        double meshSize;

        mesh_ptrtype mesh_fluid, mesh_reference_fluid;
        MeshMover<mesh_type> M_mesh_mover;

        mesh_structure_ptrtype mesh_struct;

        space_fluid_ptrtype Xh;
        element_fluid_ptrtype U_fluid;

        element_fluid_velocity_ptrtype M_identity_ale;
        element_fluid_velocity_ptrtype M_mesh_velocity;

        bdf_ptrtype M_bdf_fluid;

        structure_space_ptrtype Sh;
        element_structure_ptrtype M_depl_top;
        element_structure_ptrtype M_depl_bottom;

        element_structure_ptrtype M_bc_reference_top;
        element_structure_ptrtype M_bc_reference_bottom;

        std::vector<element_structure_type> M_referencePolyBoundarySet;

        boost::shared_ptr< ALE< convex_type > > M_aleFactory;
        ale_map_element_ptrtype M_aux;

        bdf_ale_map_ptrtype M_bdf_ale_map;
        bdf_ale_velocity_ptrtype M_bdf_ale_identity;

        uint M_cpt_export;
        double rho;
        double mu;
        double penalbc;
        double Vitesse_Entree;
        double Pression_Sortie;

        boost::shared_ptr<export_type> exporter;



    }; // Stokes
} // Life
#endif /* __NSSTOKES_H */
